A Consistent Framework for Coupling Basal Friction With Subglacial Hydrology on Hard‐Bedded Glaciers

Abstract Below hard‐bedded glaciers, both basal friction and distributed subglacial drainage are thought to be controlled by a network of cavities. Previous coupled hydro‐mechanical models, however, describe cavity‐driven friction and hydraulic transmissivity independently, resulting in a physically inconsistent cavity evolution between the two components of the models. Here, we overcome this issue by describing the hydro‐mechanical system using a common cavity‐evolution description, that governs both transient friction and hydraulic transmissivity. We show that our coupling approach is superior to previous formulations in explaining a unique observation record of glacier sliding speed from the French Alps. We find that, at multi‐day to multi‐decadal timescales, sliding speed can be expressed as a direct function of basal shear stress and water discharge, without accounting for water pressure, which simply adjusts to maintain the cavitation ratio needed to accommodate the water supply.


Introduction
In this supporting information we provide further details on the model equations and set up (Text S1), on the estimation of parameters (Text S2) and on the comparison of our model predictions with previous approaches (Text S3 and S4). This supporting information also includes additional figures and tables that are referenced in the main text.
Text S1. Model description

Ice flow model
The model solves the Stokes equation for temperate ice adopting Glen's flow law for viscous isotropic ice: where is the Cauchy stress tensor (MPa), g is the gravity acceleration (m s -2 ) and is the ice density (kg m -3 ). Basal boundary conditions are determined by the friction law described in Eq.
(1) and Eq. (7) of the main text while a stress-free condition is assumed at the surface.

Subglacial hydrology model
The subglacial hydrology model mainly uses the implementation of the GlaDS model (Werder et al., 2013) in Elmer/Ice (Gagliardini & Werder, 2018) but integrates the cavitation state to estimate the characteristics of the distributed drainage system at any given time. Distributed (linked cavities) and channelized water flows are coupled together in order to compute the hydraulic potential Φ (MPa) at the glacier bed defined as Φ = Φ + with Φ = the elevation potential at the glacier bed and the water pressure. Flow within the distributed drainage network follows the following governing equations: where Q s is the water sheet discharge (m 2 s -1 ), k s 0 is the sheet conductivity, (α s , β s ) are constant exponents, m input is the water input (m s -1 ), e v is the englacial void ratio, ρ w is the water density (kg m -3 ) and (p 1 , p 2 ) are constant exponents linking the cavitation state to the sheet thickness ℎ as ℎ = ℎ and to the hydraulic conductivity as = . Flow within channels is solved on the element edges and remain the same as described in Gagliardini & Werder (2018). Φ is set to be equal to Φ at the glacier front, and a no input flux condition is applied at the glacier outline boundary.

Coupling ice flow and subglacial hydrology
Following a previous implementation of GlaDS in Elmer/Ice (Gagliardini & Werder, 2018), the hydrological model is solved within the finite element code Elmer/Ice on a 2D mesh at the glacier base, while the Stokes equation is solved in 3D on the whole mesh. At any given time, the hydrological problem is first solved for a given sliding velocity in order to provide a solution for the hydraulic potential Φ and cavitation state . The cavitation state is then used to solve the Stokes equation and to update the sliding velocity. In order to limit the total computing time, we consider a smaller time step for solving hydrology than that for solving ice flow, given that hydraulic potential variations are much faster than the cavitation state evolution which drives the change in ice flow. We therefore use a time-step of 0.8 hour for the hydrology and 2 hours for the ice dynamics. By construction, the cavitation state results from an averaged state of effective pressure at the timescale for which the cavities evolve (several hours). There is thus no need to average the hydrological condition when driving the ice flow model at a larger time-steps as we do. We also find out that using a bigger time-step of up to 5 days for the ice dynamics does not affect the results, since the evolution of θ is not sensitive to short term velocity variations but mainly to multi-day averaged velocities. In order to further limit the computation time when doing for parameters estimation (Text S3), there we used a 100-hours' time-step for the ice dynamics.
4. External forcings: using a simplified glacier geometry and water flux description As described in the main text, we adopt a simplified geometry represented by an 0.6 X 5 km ice slab over a bedrock with uniform slope (6° slope, Figure S1). The surface topography is constructed using the "plastic approximation" which consists in ensuring uniform driving stress and thus constructing the surface elevation by solving (Cuffey & Paterson, 2010): where H is the glacier thickness (m), the surface slope and the bed slope. Eq. (s5) can be rewritten in the form: which allows to compute numerically the elevation of the surface starting from = 10 m at the glacier front ( Figure S1). The temporal evolution of the surface is then imposed at daily time scale by prescribing the temporal evolution of basal shear stress reconstructed in Gimbert et al. (2021) at the location of the sliding speed measurements on the Argentière Glacier.
The coupled problem is then solved on the given evolving geometry and is forced by a time-dependent water source term (see Eq. (s2)) prescribed at the glacier base as a function of surface elevation. We assume the source term to be equal to the surface melting (m w. eq. s -1 ) which is computed with a degree day model forced by air temperature time series and calibrated in order to have a modeled water outlet discharge being in agreement with that observed. We use daily maximum and minimum air temperatures measured at the Lyon-Bron meteorological station applying a 6.5 K km -1 lapse rate to estimate the Argentière Glacier air temperatures (Vincent, 2002). Sub-daily temperatures are computed by a sinusoidal function scaled by maximum and minimum air temperature.

Text S2. Model parameter estimation
The model parameters are summarized in Table S1 with the tuned parameters being highlighted in red. α s , β s, α c , β c , e v are taken from the literature (Gagliardini & Werder, 2018;Werder et al., 2013) and C, h r , l c are fixed to realistic values. The other parameters (A s , m, q, k s 0 , p 1 , p 2 , k c , A c ) have been calibrated to minimize the mis-fit between the model and the sliding speed record of Argentière Glacier from 2000 to 2018. The model outputs are extracted and averaged along a cross section located at 1500 m from the glacier front, which is not affected by the inflow boundary condition when solving the Stokes equation. The likelihood P of the tuned parameters is computed from the root mean square error (RMSE) between the modeled and measured sliding velocities during the period 2000-2018 assuming Gaussian distribution: The standard deviation on the RMSE is fixed to 0.3 cm d -1 .
In winter, when no surface melting occurs, we assume that glacier sliding follows the Weertman sliding law with: = (s10) where is the Weertman sliding parameter (m s -1 Pa -m ) and the Weertman sliding exponent. We first constrain and using sliding velocity and basal shear stress as measured in winter when sliding velocity is minimum (using a one-month averaged value centered on the minimum seasonal velocity each year). We find the best accordance with data for = 6.8083·10 -26 [5.6 -8.4] ·10 -26 m s -1 Pa -3.72 and = 3.72 ± 0.24 using the MATLAB ® Curve Fitting toolbox ( Figure S7). In the following, these values are fixed and prescribed in the coupled model prior to calibrate the other parameters.
We find that the summer speed-up follows the theoretical relationship demonstrated in this study between discharge and cavitation state, which suggests the channel system plays a minor role during the melting period at the Argentière Glacier (Figure 2). The calibration of the channel parameters ( , ) appears to be relatively independent of that of the cavity network ( , , , ). We therefore first find the best parameters for the channels network with fixed values of ( , , , ) that gives reasonable agreement between model and data in the early melting season. Using the whole sliding velocity time series (2000-2018), we find a best match between model predictions and observations for A c = 6.0 · 10 -24 Pa -3 s and k c = 3.4 · 10 -4 m 3/2 kg -1/2 ( Figure S8). The inferred value of k c is several orders of magnitude lower than that proposed in the literature (Chen et al., 2018;Hewitt et al., 2012;Werder et al., 2013), which may be an artefact of the simplified geometry we use or due to an inaccurate formulation of the channel opening term in R-channel theories. For instance, the hypothesis that energy potential losses are purely converted in channel melting (Werder et al., 2013) may not be apply for grounded channels for which energy loss can also occur at the water-rock bed interface and thus not contribute to melt. It appears that the best values of (A c , k c ) are the ones that do not allow channels to develop significantly during the melting period but still allow channels to reduce the water pressure in winter when the sheet conductivity vanishes. The development of significant channels in summer systematically would lead to low sliding velocities in August/September that are in disagreement with the data. However, channels are still needed to explain the sliding velocity in winter as they allow to decrease water pressure and sliding speed during this period in better agreement with the data. In winter, channels exist as a dense network of small conduits and can be viewed as melt-driven opened conduits that maintain connections between cavities when mechanical, cavitation-driven opening is not sufficient.
For fixed values of A c and k c , we then explore the parameter space ( , , , ) by accounting for the fact that the values of and are constrained by the observed relationship between discharge and sliding speed, which gives + = 7.2 (Figure 2a). We thus only optimize the value of (p 1 , k s 0 ) for different values of . We find that = 1.5 best explains the data ( Figure S2a) and that the value of q is not affecting the best values of the couple (p 1 , k s 0 ). We obtain p 1 =0.6 and k s 0 = 0.1 m 7/4 kg -1/2 (Figure 2a) leading to =6.5. The low value of p 1 and the high value of p 2 reflect the fact that cavities open quite fast but remain unconnected until a certain threshold close to the limit of stability (transition strengthening/weakening) is reached ( Figure S2b). We note that the value of C does not affect the results and only has the consequence of shifting the effective pressure in a way that the quantity C*N is conserved.

Text S3. Details on Eq. (13) of the main manuscript
We combine equations ( Gagliardini et al. (2007) friction law gives (see Eq. (2) of the manuscript): and combined with Eq. (s8), it gives: Introducing the expression of h from Eq. (s7) gives Eq. (13) of the manuscript:

Text S4. Comparison with the original Elmer/Ice -GlaDS coupled model (Gagliardini and Werder, 2018)
To compare the performance of the present framework with the previously proposed ones in explaining our data, we perform simulations using the original coupling between GlaDS and Elmer/Ice (Gagliardini & Werder, 2018). In this simulation, we however keep the transient friction approach as introduced in Thøgersen et al. (2019) (with d c =0.5 m) in order to be able to use rate weakening friction and hourly water input. Conventional steady friction laws would indeed prevent the model to converge as effective pressure often exceed the threshold of rate weakening friction ( / ) when providing water input rate at hourly time step. The coupling is done in a standard way with GlaDS solving for effective pressure at a given sliding velocity and independently of the frictional state. The sliding velocity is then updated by solving the ice flow problem for a new given effective pressure. This simulation allows to test the influence of imposing the consistency between friction and hydraulic transmissivity as done in the present approach. We optimized the distributed system parameters ( , and ) to obtain the best fit with the sliding velocity observations. The other parameters are set to the one used in our study. Results are described in the main manuscript with reference to Figure 2 and Figure S4. Figure S1. Idealized geometry used to solve the coupled problem. The dimensions are 600m width and 5km long for a mean bed slope of 6 degrees in order to approximate the Argentière Glacier setting. Figure S2. Parameter likelihoods computed from the misfit between modelled and measured sliding speed using the consistent hydro-mechanical coupling introduced in this study. The panels in (a) explore the parameter space (p1, k s 0 ) for different values of q and with other parameters fixed at the reference values given in Table S1. The panel (b) shows the function h(θ) and k s (θ) obtained for the best fit parameters with q=1.5.  Figure S4. Root mean square error between modeled and observed sliding velocity during the period 2000-2018 using a standard coupling approach between Elmer/Ice and GlaDS (Gagliardini & Werder, 2018) for different values of the flow rate factor (closing term of cavities) compared to the error obtained using the consistent hydro-mechanical coupling (upper left panel).

Figure S7.
Observed sliding velocity at the Argentière wheel (red) and modelled sliding velocity using a Weertman law calibrated on observed minimum velocities each year. See section S2 for basal shear stress estimation used here. Figure S8. Channel parameters likelihood computed from the misfit between modelled and measured sliding speed using the consistent hydro-mechanical coupling introduced in this study. The figure explores the parameter space (A c ,k c ) for other parameters fixed at the reference values given in Table S1.  Channel discharge Q m 3 s -1